Low vorticity and small gas expansion in 
premixed flames 

Bruno Denet^ and Vitaly Bychkov^ 
9th February 2008 

^IRPHE 49 rue Joliot Curie BP 146 Technopole de Chateau Gombert 13384 
Marseille Cedex 13 France 

^Institute of Physics, Umea University, SE-901 87, Umea, Sweden 

submitted to Combustion Science and Technology 
Abstract 

Different approaches to the nonlinear dynamics of premixed flames 
exist in the literature: equations based on developments in a gas ex- 
pansion parameter, weak nonlinearity approximation, potential model 
equation in a coordinate- free form. However the relation between these 
different equations is often unclear. Starting here with the low vor- 
ticity approximation proposed recently by one of the authors, we are 
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able to recover from this formulation the dynamical equations usually 
obtained at the lowest orders in gas expansion for plane on average 
flames, as well as obtain a new second order coordinate-free equation 
extending the potential flow model known as the Frankel equation. It 
is also common to modify gas expansion theories into phenomelogical 
equations, which agree quantitatively better with numerical simula- 
tions. We discuss here what are the restrictions imposed by the gas 
expansion development results on this process. 

1 Introduction 

The nonlinear description of the Darrieus-Landau (DL) instability of pre- 
mixed flames began twenty-five years ago, when Sivashinsky obtained, as a 
first-order development in powers of a gas expansion parameter, a nonlinear 
equation known today as the Sivashinsky equation (Sivashinsky 1977). Start- 
ing with the first simulations of Michelson (see for instance (Michelson and 
Sivashinsky 1982)) and with the analytical solution of Thual with coauthors 
(Thual, Frisch and Henon 1985), this equation has shown a surprising quali- 
tative agreement with experiments and direct numerical simulations (Denet 
and Haldenwang 1995, Bychkov, Golberg, Liberman and Eriksson 1996, 
Kadowaki 1999, Travnikov, Bychkov and Liberman 2000). The only nonlin- 
ear term of the Sivashinsky equation has a purely geometrical origin and is not 
related to the usual nonlinearity of the Navier-Stokes equations. Challenged 
by Clavin about the fact that the Navier-Stokes nonlinearity should induce 
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modifications of the flame equation, particularly for realistic gas expansion, 
Sivashinsky went on to show, in a joint paper with Clavin (Sivashinsky and 
Clavin 1987), that even at the second order in gas expansion, the equa- 
tion obtained was still an equation with the same terms, only with modified 
coefficients. Today various methods try to improve on these equations by 
using the approximation of weak nonlinearity or next orders in gas expan- 
sion (Zhdanov and Trubnikov 1989, Joulin 1991, Bychkov 1998a, Kazakov 
and Liberman 2002b, Boury 2003). On the other hand, similar to the origi- 
nal Sivashinsky equation, all these equations have been derived for the pla- 
nar on average flame front, and are valid only when the slope is not too 
large. Actually, different variations of the Sivashinsky equation have been 
constructed for different geometries like expanding flames (Filyand, Sivashin- 
sky and Frankel 1994, D'Angelo, Joulin and Boury 2000) or oblique flames 
(Boury and Joulin 2002), when the flame shape departs slightly from the 
unperturbed case. 

A different equation was introduced in the theoretical community by 
Frankel (Frankel 1990) in 1990, although parts of this approach were antici- 
pated by numerical studies, see for instance (Ghoniem, Chorin and Oppenheim 
1982). The Frankel equation was derived in a coordinate-free form (in two 
or three dimensions) using similarities between electrostatics and the flame- 
generated flow. Constructing his theory Frankel assumed that the flow is 
potential everywhere and neglected vorticity generated by the flame in the 
burnt gas. As a matter of fact, such an assumption originated in the analy- 
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sis by Sivashinsky (Sivashinsky 1977). Numerical simulations demonstrated 
that the Frankel equation describes qualitatively well the nonlinear evo- 
lution of expanding flames (Frankel and Sivashinsky 1995, Blinnikov and 
Sasorov 1996, Ashurst 1997, Denet 1997) and oblique flames (Denet 2002) 
(even recently in three dimensions (Denet 2004)). The Frankel equation be- 
came especially popular in the studies of fractal flames developing because of 
the DL instability on large scales (Blinnikov and Sasorov 1996, Denet 1997). 

From a theoretical point of view, the relationship between the Frankel 
and Sivashinsky equations was put forward from the very beginning, since the 
original paper (Frankel 1990) showed that the Frankel equation reduces to the 
Sivashinsky equation for plane on average flames (with lateral boundaries at 
infinity) and for expanding fiames. However, validity of the Frankel equation 
has been questioned rather often because of the assumption of the potential 
flow, which, in principle, violates hydrodynamic conservation laws of the 
flame front. The original analysis (Frankel 1990) did not demonstrate if the 
Frankel equation follows from the hydrodynamic equations in the limit of 
small gas expansion, as it was done for the Sivashinsky equation (Sivashinsky 
1977). 

Recently, however, one of us (Bychkov, Zaytsev and Akkerman 2003) re- 
considered the problem and introduced a low vorticity approximation (com- 
pared to the Frankel case, where the vorticity is strictly zero), which was 
justifled by the previous analysis of curved flames (Bychkov 1998a). The 
approximation of low vorticity enabled to derive a system of coordinate-free 
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equations describing the evolution of the front even for realistically large 
thermal expansion of the burning matter. This system of equations is rather 
complex and has not been successfully solved numerically for the moment. 
In the present paper, we develop this system up to the second order in gas 
expansion, which is equivalent to developing the hydrodynamical equations 
at this order. Such calculation proves the Frankel equation at the first order 
in gas expansion, and leads to a second order form of the equation (we will 
however insist on some difficulties specific to the oblique fiame geometry). 
This form in turn can be demonstrated to contain the Sivashinsky-Clavin 
equation in the planar on average case. 

In Section El we obtain this second order form. In Section [HI this equation 
will be reduced to the Sivashinsky-Clavin equation. In Section |3| we discuss 
some basic problems related to the Sivashinsky-Clavin equation and expan- 
sion in powers of a small parameter in general. Finally, Section El contains a 
conclusion. 

2 Derivation of the second order Frankel equa- 
tion 

Let us start this section by summarizing the main points of the low vorticity 
approximation of (Bychkov et al. 2003) (the reader is invited to read this 
article for more details). This approximation is derived from the hydrody- 
namical equations : 
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V-u = 0, 



— + (u • V) u = -i?vn, 

OT 

where the equations written are made non-dimensional with the use of the 
laminar flame velocity Uf and a reference length R. The non-dimensional 
velocity is noted u and the pressure H = [P — Pf) /pj-Uj, — 1 in the 
fresh mixtures and {} = Q in the burnt gases. The flame is considered as a 
discontinuity, © = Pf/Pb is the ratio of density in fresh and burnt gases. The 
parameter (© — 1) will be the parameter of the expansion. The boundary 
conditions of the flame can be classically shown to lead to 

u+ = u_ + (© - 1) n, 

n+ = n_ + 1 - ©, 

where n is the normal vector to the flame surface, directed towards the burnt 
gases. The first boundary condition accounts for both the jump of normal 
velocity and conservation of tangential velocity at the front. We introduce 
the velocity potentials in fresh (— ) and burnt (+) matter ^-t, which satisfy 
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the Laplace equation and are defined by 

u_ = V<^_, 

u+ = Up+ + u„+ = V</)+ + u„+, (1) 

where Up+ and Ui;+ are the potential and vortical parts of the velocity field 
in the burnt gases. 

We work in a local system of coordinates moving with the fiame front at 
the velocity — V^^n such as Vs = ej ■ V, where is the unit tangential vector, 
and 

d d ^ d 
Ots dr ^ dn 

Using basic properties of Green functions of the Laplace equation, Bernoulli 
integrals, and boundary conditions at the front, the following system of equa- 
tions (low vorticity limit) is obtained 

— (0+ - 0_ + 0„ (1 - 6)) = l^u^_-{Q - 1) v;2+(0 - 1) Vs+VsU^r.++f, 

(2) 

^ + (up+ ■ V) = 0. (3) 

Note that the form given here is not the final form of (Bychkov et al. 2003), 
here we do not incorporate G into the variables in order to make the devel- 
opment in (B — 1) easier. Equation ^ comes from the Bernoulli integral; 
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/ is generally a function of time appearing in the Bernoulli integrals (unlike 
(Bychkov et al. 2003) we have included a constant term (9 — 1)^/2 into /). 
We will choose however not to include any additive terms containing a func- 
tion of time in the potentials, so that / has to be considered as a constant. 
Furthermore, strictly speaking, we could add a constant in u„+ and subtract 
it from Up+ but naturally u„+ has to be taken as small as possible, which 
makes the Bernoulli integral a good approximation of the Navier-Stokes equa- 
tion in the burnt gases. A non-zero value of / would lead to a constant value 
of Uvn+ at infinity (convected by the potential velocity) for a plane or oblique 
geometry, so that we must have in this case / = 0. Similarly, in the case of 
the expanding geometry, a non zero value of / would leave a constant value 
of Uyn+ behind the front, which is not possible because there is no source or 
sink present inside the expanding fiame. 

Equation ^ expresses the fact that the vortical part of the velocity field 
is convected by the potential part (it will create a shear fiow at infinity in 
the plane configuration that would disappear far from the fiame if viscosity 
was included). Note that this equation is different from the equivalent one 
given in (Bychkov et al. 2003), which included only the time derivative. As 
formulated in (Bychkov et al. 2003), the approach of low vorticity takes into 
account only convection by the uniform component of the velocity field. In 
the geometry of planar (on average) fiame front propagating along z-axis 
this corresponds to the drift term Qduy/dz. In the case of expanding fiames 
there is no uniform velocity in the burnt matter and the vorticity created at 
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the flame surface is simply left behind as the flame radius increases; by this 
reason the drift term was omitted in (Bychkov et al. 2003). In the present 
paper we are interested in small gas expansion of the flame equations, which 
allows to consider a more general form of equation ^ . Note also that at the 
dominating order in gas expansion, the convection by the potential velocity 
is simply equivalent to the convection by the injection velocity. Developing 
equation ^ and using the continuity equation (n is the normal, s represents 
the tangential coordinates) 



— h Vs ■ u„j+ = 

on 



and u^t+ = Vs (0_ — 0+), we obtain 



n + {{upt+ ■ Vs) ^v) ■n = {Vs + Upn+) (0_ - 0+) , (4) 



where K = 1 — Compared to (Bychkov et al. 2003), we have two 

supplementary terms coming from the convection by the potential flow field. 
The first term of equation ^ (the term with the time derivative) is also 
slightly corrected compared to this paper. 

The Laplace equation leads to an equation with a different form in two 
and three dimensions: 



3D : = f( ® + (0+ - 0-) n . dS{r, 



2n 



r 



(5) 



or 



2D : <P++4>- = -- / ( (0 - 1 - Uyn+) In Ir^ - r| - (0+ - 0_) n ■ ^ ] dlivs 

TT J \ r« — r / 



(6) 

Now let us expand all variables in powers of (6 — 1): 
0± = + + • • • 

(1) , (2) , 

V, = + vp + vp^ ■ ■ ■ 

_ (1) (2) 

^ = 0(6-1), 

where the subscript means that the term is of order (6 — 1)*. Note 
also that we do not expand the positions r in powers of (6 — 1); this can be 
done only in particular geometries, when the difference between the actual 
and the unperturbed positions is of order O (6 — 1). To reduce the equation 
obtained to the Sivashinsky-Clavin equation in the planar geometry, we will 
have to consider the fact (see sectionlHl) that the vertical coordinate is of order 
(6 — 1), but this order of magnitude estimate is not geometry-independent. 

First, using the relation K = 1 — we find that = 1. Note that at 
this zeroth order, we can have a term ^foundary = ^inj (a constant injection 
velocity, for instance) that has to be added to the velocity field in order to 
satisfy the boundary conditions at infinity. At each order, we will encounter 
a ^boundary term, SO let us define it properly, ^boundary is a velocity field 
obeying the Laplace equation, without jumps on the fiame surface, which 
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helps satisfying the boundary conditions at infinity. Although we will perform 
the calculations in a reference frame without injection velocity, let us note 
that with an injection velocity we would have = 1 — Vinjn- Obviously, 
the Bernoulli relation is Galilean invariant, so we choose the reference frame 
where the calculations are simpler, knowing that at the end, the injection 
velocity (if present) may be added to the final formula. 
By developing equation Q we have at the 0(1) order 

- (6 - 1) + (6 - 1) + V^'^u^l = 0, 

which, with Vj^'^ = 1, leads to u^n+ = 0. We also have, according to equation 

(H 

= (0« - 0«) , 
which implies in three dimensions, using equation ^ and the value of u^n+ 

This is the potential of a uniformly charged surface. Apart from multiplica- 
tive factors we will call (6 — 1) the surface charge (or line charge for the cor- 
responding two dimensional equation). The reader is referred to (Denet 2002) 
for a simple presentation of this electrostatic analogy. This potential leads 
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to the corresponding velocity 

Vm = = (l + i- / !i-|£lZi:)<iS(r.)) - V«„^,, . (8) 

Up to the first order, Vg = 1 + Vj^\ which is exactly the Frankel equation in 
three dimensions (when the flame front is a surface). 

Now let us consider the 0(2) terms. At this order, we have from equation 

m 

-2 (9 - 1) Vj'^Vj'^ + (6 - 1) Vj'^ + + K^^^^l = 0. 

With the previous computed values V^'-^-' = 1 and = 0, it reduces to 



u 



(2) 
vn+ 



(e-i)i^«. (9) 



Taking into account that 0^"* = 0+"* equation (jH) gives 

n, (10) 



V2 /a(2) _ a(2) 



^U-mj ■ Vs) ^ 



(2) 



SO we obtain from ^ in three dimensions with the previous value of u. 



with (0^'' — 0!^'') determined by equation (fTn|) . Let us note however that 



12 



this term exists only when there is a strong tangential velocity field (i.e. 
for oblique flames). In the plane and expanding flame cases, the tangential 
velocity field can be neglected, so that [(t)L^ — 0+"*) is also negligible in the 
previous formula. In the rest of the article, we will only write the formulas 
with this term neglected, but let us insist on the fact that by doing so, we 
neglect a dipolar contribution to the potential, which could be important in 
oblique cases. 

Summing the first and second order terms, we have 

+ ^ A" ^'^f-rJ (") 

which can be interpreted as a local surface charge modified from © — 1 to 

©-i-(e-i)yw. (12) 



Then the total velocity field up to the second order is 



,(1) .,(2) 



2 1^ 27ry |r, - 



|3 



dS{vs 



r| 

boundary boundary boundary i \ ) 
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or using V^^^^ = K — 1 + 0(0 — 1) we find the equivalent formulation 
Vs = l + ^^(2 -Vs)(l + ^ I ^^-^^m^rfS(r,)') - \,ound^ry " n. (14) 



Formula (fTTS|l is the coordinate-free equation for the flame front velocity ob- 
tained within the second order in gas expansion; the first order gives the 
Frankel equation. We recall that Vg is the normal velocity of the front prop- 
agation, including laminar fiame speed and induced velocity field. Curvature 
and strain effects are not included, see (Bychkov et al. 2003). In two dimen- 
sions, the equation, obtained by similar calculations is 



VT J r, — r 



— V^°^ . n — V^^^ ■ n — V*-^^ ■ n 

* boundary boundary boundary i v-"-"/ 



or 



— 1 / Irn-fr— \ 

Vs = l + ^^(2 - K) 1 + - / , ^ ' J dl{rs) - Y,,.,r,dary " H. (16) 

In order to solve this equation numerically, the front would have to be de- 
scribed with marker particles moving according to the equation 
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Numerical solution to equation (fT7j) may involve reconnections in two and 
three dimensions; possible ways to overcome these difficulties are described 
in (Denet 2002, Denet 2004). 

3 Reduction to the Sivashinsky-Clavin equa- 
tion 

We have previously said that equation (fT5|l obtained in the coordinate-free 
case is supposed to be the equivalent of the Sivashinsky-Clavin equation in 
the plane case. Indeed, this last equation is derived by a development which 
is one order higher than the Sivashinsky equation, while equation (|T5|l is one 
order higher than the Frankel equation. However, transition from one case 
to the other is not obvious at all. The Sivashinsky-Clavin equation contains 
the same terms with different coefficients, but our modified Frankel equation 
contains a surface (or line) charge, which is not apparently constant at every 
position on the front. Furthermore, the order of magnitude of the flame 
front position makes another problem, since the position is supposed to be of 
order O (0 — 1) for the planar case and 0(1) in the coordinate-free case. So, 
is it possible to recover the Sivashinsky-Clavin equation from equation (fTKll ? 
Naturally, we have the help of the original Frankel article (Frankel 1990), 
where it was shown for the planar (on average) case with lateral boundaries 
at infinity and for the circular expanding case that the Frankel equation 
reduces to the Sivashinsky equation. Let us start by repeating this reasoning 
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before going to the next order. In two dimensions, the Frankel equation gives: 



•^s — -L -I- 2 y-^^J Ij. _ j.|2 "H^^^sJ I ^boundary " ^boundary 



In the planar case, the normal vector is 



n(r) = [-ay,l]/^l + al, 



where y is the lateral coordinate, a is the vertical position of the front, the 
vertical z coordinate is positive towards the burnt gases, ay is a notation for 
the y derivative. Then equation (fT7|l gives 



-Vs. (18) 



1 + al 



In order to recover the Sivashinsky equation, we have to develop the square 
root, so we suppose a = O (6 — 1) and y = 0(1). We have 



al 9- 1 e-1 /■+°° a(x + y) -a{y) - ay {y)x , yio) yW 

boundary ' ^ z boundary ' 



Frankel has shown that the principal value in this formula is 



where I (a) is the Landau operator (multiplication by \k\ in Fourier space. 
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see (Frankel 1990) for details). 

Let us take as usual the hypothesis that the injection velocity is parallel 
to the z direction and has the value ^f^undary ~ ^ inj — [0, 1] . This term 
simplifies itself with the —1 of the previous equation. All the following terms 
of the boundary velocity are obtained by saying that they do not change the 
injection velocity. If this velocity is imposed at a finite distance, then this 
is simply obtained by the same integral as in the Frankel equation, with the 
same charge, but for the mirror image of the front relative to the injection 
location (see (Denet 2002) for an example). If the injection is moved to 
infinity, then the velocity field, according to the Gauss theorem, will be the 
same as the velocity generated by a plane front, but with a charge multiplied 
by the ratio of the front surface to the surface of an equivalent plane front 
y'^^fyandary — + "^y) — ^f^' whcrc () is a lateral mean value. Higher 

orders will be obtained by developing the square root and replacing — 1 by 
— 1 — (0 — 1) V^*^^\ as discussed in the previous section. For the time being 
however, the boundary term only leads to suppress the term — (0 — 1) /2 . 
We thus obtain the Sivashinsky equation (without curvature terms) 

+ f = («) ■ (19) 

Note that with a = O (0 - 1) and g = O (0 - 1) all the terms of the 
equation are of the second order. 

Let us now perform the same calculation at the next order starting from 
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the two-dimensional equation (fTK|l 

^ ^ ^ e-i-(e-i)K'" A ^ W ILi£^,,(„) . 

2 \^ ttJ |rs-r| / 

Equation (fH^ is still valid. We have, retaining terms up to the third order 

_ _,_gl_9 - 1 - (e - 1) V.'" _e-2^.^9 - 1 - (9 - 1) vy 

2 2 4 ^ 2 ^ ^ 

As before, the constant terms (not depending on y) are eliminated by intro- 
ducing the boundary velocity terms, but what do we do about V}^\ given in 
equation (jH]), which depends on position? The answer is simple: since now 
we have a = O (0 — 1), then we obtain a product of — 1 and a in V"^*^^-*. 
This term is not of the first order anymore, it is now of the second order and 
can be neglected (note that this is not necessarily the case in all geometries). 
The constant terms in Vj-^^ are as usual suppressed by the boundary velocity, 
and we find 

al 0-l-(0-l)((0-l)^) 0-1 ^ 0-l-(0-l)((0 
= -^-f 2 ^"^+ 2 

Finally, adding the effects of the boundary terms, and disregarding terms of 
higher orders we obtain 

al f 0-l\ &~l/al\ 0-1/ 0-i\ 
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This is the Sivashinsky-Clavin equation, written in units of laminar flame 
speed (relative to the fresh gases). The lateral mean value term comes from 
the boundary velocity and was absent in the original article. It was later 
added by Joulin, who used it in a number of papers, the first one being prob- 
ably (Joulin and Cambray 1992) and called it a counter-term. In Sivashinsky- 
Clavin units the value 7 = was used to characterize gas expansion and 
all velocities were scaled by the laminar flame speed relative to burnt gases 
Ub = QUf = t///(l — 7). In such units equation (f2n|) would be written as 




In units of laminar flame speed relative to fresh gases Uf, but using 7 as a 
parameter the equation takes the form 

4 Do we have the first or second time-derivative 
in the flame equations? 

One more question concerns the order of time-derivative in the Sivashinsky- 
Clavin equation, which also affects the structure of nonlinear terms in the 
equation. The original DL dispersion relation (Landau and Lifshitz 1989) 
is of the second order describing two independent linear modes of the flame 
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front perturbation: one mode is growing and one is decaying. Unlike this, 
the Sivashinsky equation is of the first order in time capturing the growing 
mode only; the Sivashinsky-Clavin equation has the same property. When in- 
vestigating development of the DL instability, the growing mode dominates 
over the decaying one and the first-order derivative looks quite sufficient. 
However, the second-order time derivative proved to be of principal impor- 
tance in many adjacent problems of flame dynamics: propagation of tulip 
flames (Dold and Joulin 1995), flame interaction with sound (Searby and 
Rochwerger 1991, Bychkov 1999), with shocks (Bychkov 1998b) and with 
external turbulence (Searby and Clavin 1986, Aldredge and Williams 1991, 
Akkerman and Bychkov 2003). In all these cases some external force redis- 
tributes energy between growing and decaying modes, and it is incorrect to 
exclude the decaying mode out of consideration. Besides, as we said above, 
changing the order of the time derivative we also modify the nonlinear terms 
of the equation. At this point, following comments of one of the referees, 
we would like to add that, in general, even a first-order time derivative may 
produce a complicated spectrum with many stable and unstable modes. How- 
ever, in the particular case of the linear DL instability for an infinitely thin 
planar flame front and a flxed wave number of perturbations the number of 
modes is unambiguously related to the order of time derivative. The flrst 
order time-derivative of the Sivashinsky equation provides only one mode 
(growing), while the original DL dispersion relation with the second order 
derivative has two modes (one is growing and one is decaying). 
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Let us consider how the Sivashinsky-Clavin equation should be modified 
to take into account the second time derivative. Within the accuracy of the 
Sivashinsky-Clavin approach we can substitute (fm into the second-order 
terms of (1201) and find 



e-1 62 /«?A 0-1 



at + at + -al + {^) = ^—I{a) (21) 



or 



r I, ^ r l(dOLl\ e 2 0-1/ 2\ ^'^r, ^ 



The operator has a meaning of an integral, and it is defined with the ac- 
curacy of a constant. Still, this constant is included already into the counter- 
terms. One more trouble with this operator concerns acting on a con- 
stant. However, in the case of equation (f^ . a constant under would 
imply physically meaningless solutions like a planar accelerating fiame front, 
which may be obviously ruled out. One can easily see that the linear terms 
of equation (f22|l coincide with the DL dispersion relation written in the limit 

of e - K 1 

^I-'M+a, = ^I{a). (23) 

Equation contains only one counter-term, because the time-dependent 
nonlinear term involves the complete time derivative and gives zero after 
averaging. Though equations (f2n|l and (|22|l are mathematically equivalent 
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within the expansion in 6 — 1 ^ 1, they may lead to somewhat different 
conclusions about properties of curved flames. For example, one of the most 
important questions in the nonlinear theory of the DL instability is the ve- 
locity increase of curved stationary flames. Let us consider propagation of 
such a flame a{t,y) = —Qt + a{y). In that case equations (f2n|) and (f22|l 
reduce to 

and 

-n + -al- ^- {al) = ^-/(«). (25) 

As we can see, equations and are different. It is interesting that 
equation (f25|l is consistent with the stationary theory (Bychkov 1998a) within 
the accuracy of (0 — l)a^, while ([2^ is not. The above calculations illustrate 
the fact that any rigorous expansion in power series leaves plenty of freedom 
for mathematical manipulations, which may lead to ambiguous physical con- 
clusions. Unfortunately, almost always people use expansion in power series 
of a small parameter to obtain physical results beyond the validity limits of 
the expansion. As a simple illustration, suppose that we calculated some 
value a using expansion in power series of e <^ 1, and found a = 6 in the 
zero order approximation with a = 6(1 — ce) for the first order. The same 
expression may be written in an infinite number of equivalent mathematical 
forms like a = b{l — [c + b — a]e), a = b{l — 2c£:)^/^, a = 6/(1 + ce), etc. 
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Though these forms are equivalent within the first order in £ <S 1, we come 
to different conclusions when investigating zero points of a with the help of 
these expressions. In the above four versions of the same formula we find 
a = at £ = 1/c, e — l/{c + b), e — l/2c and £ = oo, respectively. One 
encounters a similar trouble within both linear and nonlinear theories of the 
DL instability. For example, the linear theories (Pelce and Clavin 1982) and 
(Matalon and Matkowsky 1982) lead to noticeably different expressions for 
the cut off wavelength of the DL instability, though both theories are math- 
ematically correct and equivalent within the same accuracy of small wave 
numbers. Performing manipulations similar to those described above, we 
can actually obtain infinite number of absolutely different formulas for the 
cut off wavelength in scope of the theory (Pelce and Clavin 1982) keeping 
the same accuracy. So in the case of the linear DL instability the simple 
example a = b{l — ce) is sufficient to explain the discrepancy between the 
analytical values for the cut-off wavenumber obtained by different authors, 
as this formula is a simplified version of the dispersion relations obtained in 
the two previously mentioned articles. 

What does it mean if one tries to derive a non linear equation for the 
DL instability using perturbation methods ? In that CclSGj clS pointed out by 
one of the referees, the perturbations concern operators instead of functions, 
which is a much subtler subject. Of course, in this case, the example with a = 
b{l — ce) is ultimately simplified, still it gives a rough idea why the non linear 
approaches like (Zhdanov and Trubnikov 1989, Joulin 1991, Bychkov 1998a, 
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Kazakov and Liberman 2002b, Kazakov and Liberman 2002a) may lead to 
different equations for a flame front even if all mathematical calculations are 
performed correctly. Singular perturbation methods for partial differential 
operators are much more sophisticated than for functions, and can involve a 
large number of new phenomena, boundary layers being of course the most 
well-known example. However singular layers can also occur during the time 
evolution, for instance initial layers for times close to the initial conditions. 
See below for a discussion of the difficulties that could happen in formal 
manipulations of second order in time equations. 

When one uses expansion in powers of a small parameter, the nonlinear 
equation for a flame front may be presented in an infinite number of forms. 
For example, suppose that we have derived a time-dependent nonlinear equa- 
tion within the approach of weak nonlinearity as it was done in (Zhdanov and 
Trubnikov 1989, Kazakov and Liberman 2002b). Within the same accuracy 
of calculations one can take square of the DL dispersion relation (f23| with 
any coefficient and add it to the equation obtained. However, when we use 
the new version of the nonlinear equation to study curved stationary flames, 
square of the right-hand side of (f23|l makes a non-zero contribution, while 
square of the left-hand side becomes zero (the flrst term gives zero exactly, 
and the second term provides nonlinearity of the fourth-order). Making such 
manipulations one comes to a stationary equation like 



^+{1 + Ci) al + C2 [/(«)]' -Ci {al)-C, {[I (a)]') 



e-1 



/(a). (26) 
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with almost arbitrary coefficients Ci and C2. The only restriction on the 
factors Ci and C2 is that they should tend to zero sufficiently fast as 6 — 1, 
and that the development in 9 — 1 of the stationary solution is the same as 
the original It gives Ci = (6 - l)/2 + 0(0 - l)^ and C2 = 0(6 - l)^ 
for small (0 — 1), but it must be admitted that these restrictions are too 
loose taking into account realistically large values of 9 = 5 — 8 (also, it 
must be remembered that curvature-related terms have to be added to (|26|l 
in agreement with the linear theory of the DL instability). This result is 
rather discouraging, because it leaves no hope to obtain an unambiguous 
formula for the flame velocity by using perturbation theories. As an illus- 
tration of this fact, the authors of the two companion papers (Kazakov and 
Liberman 2002b) and (Kazakov and Liberman 2002a) have produced two 
different formulas for the flame velocity. The first-order approximation in 
(Zhdanov and Trubnikov 1989, Kazakov and Liberman 2002b, Boury 2003), 
contrary to the reasoning used above to obtain equation (f26|l . employed the 
DL-dispersion relation for the growing mode only 

Using ()27|) within the second-order approximation one can always convert 
time-dependent nonlinear terms into time-independent terms and vice versa. 
We would like to stress that manipulations with the time-derivatives is not 
something that we have invented in the present paper. They have been 
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performed already in a number of papers on flame dynamics, leading to am- 
biguous physical results. In the present paper, we just clarify the ambiguity 
and point out the danger of such manipulations. 

At this point one of the referees asks, if it is not safer to transform second- 
order time derivatives into space-derivatives. Unfortunately, even in that case 
some ambiguity remains in the non linear terms. Besides, as we pointed out 
above, by getting rid of the second-order derivatives one loses the possibility 
to study a large number of effects like tulip flames, flame interaction with 
shocks and many other phenomena. At present, we do not know how to avoid 
the ambiguity in the perturbation theory of the non linear equation for a 
flame front. In the present section we rather formulate a question than give an 
answer. We hope to insist here on the fact that an infinite number of different 
non linear equations for a flame front can be obtained by high order formal 
perturbation methods. Among these formally equivalent equations, some 
will give bad quantitative results (particularly if the development parameter 
is not small, which is unfortunately the case for flame fronts). Some will 
even give bad qualitative results, for instance we could have second order 
in time equations which, without forcing, do not approach a first order in 
time dynamics (which is well known to attract the dynamics for a flame 
without forcing). Even worse, we could have equations with pathological 
mathematical properties (this is particularly possible with time derivatives 
in the non linear terms). Although in the rest of the paper, we insist on 
the quantitative agreement on the flame velocity, it must be kept in mind 
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that at some point, the time evolution of the proposed models must also be 
compared with direct numerical simulations and found satisfactory. 

On the other hand, without using perturbation methods one comes to 
a rather complex set of equations (Bychkov et al. 2003) with almost zero 
hope to solve it analytically, and also very difficult to solve numerically. 
Luckily, in the case of curved stationary flames the problem of flame velocity 
has been solved with the help of direct numerical simulations (Bychkov et 
al. 1996); later calculations (Kadowaki 1999, Travnikov et al. 2000) conflrmed 
the original results. 

The uncertainty in the rigorous perturbation nonlinear theories of the 
DL instability increases the role of simple phenomenological models like that 
proposed in (Joulin and Cambray 1992). Using the models one can obtain 
qualitative or even semi-quantitative understanding of flame dynamics, which 
may be checked and corrected quantitatively in direct numerical simulations. 
The model (Joulin and Cambray 1992) included flrst-order time-derivative 
similar to the Sivashinsky equation (jTHI). When second-order derivatives are 
important, similar model can be constructed on the basis of equation (f22|l . 
Comparing equations (f22|l and (f23|) we can easily extrapolate (f22|l to the case 
of realistic 9 as 

^'-^ ("..) + + '-^ (If) + il + c.) - c, ^ ^/(«). 

(28) 

Similar equation was proposed in Boury's thesis (Boury 2003), in the spirit 
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of the phenomenological theory used in (Dold and Joulin 1995) to study tulip 
flames. The linear part of (f28|l is a well-known DL dispersion relation. Let us 
note however that including curvature effects in this type of equation may be 
non trivial, because actually the Markstein lengths are frequency-dependent 
(see (Joulin 1994, Clavin and Joulin 1997, Denet and Toma 1995)). The 
unknown coefficient Ci may be adjusted by using direct numerical simulations 
of curved stationary flames. In (Boury 2003) the coefficient was chosen to 
provide the same stationary amplitude as either a third order gas expansion 
theory or direct numerical simulations. However, we believe that fitting 
direct numerical simulations is a better idea. As explained before, a high 
order perturbation theory can be written in several equivalent ways, with 
different quantitative results for large 6 (see for instance, equations ([2^ and 
(f25| V Below, we illustrate that this kind of fit can actually be used to obtain 
almost any curved fiame velocity, with only some restrictions for small gas 
expansion. The stationary version of (f28|) is 

- + + Ci) - C, (al) = ^/(«). (29) 

Equation (f29|) may be solved analytically (Thual et al. 1985, Joulin and 
Cambray 1992), which leads to the maximal velocity increase 

. (^-^)\. (30) 

8(l + 2Ci)' ^ ^ 

The maximal velocity increase obtained in direct numerical simulations (Bychkov 
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et al. 1996, Kadowaki 1999, Travnikov et al. 2000) is plotted in figure [T] 
by markers. The curves of figure Q] show the analytical formulas for the 
velocity increase, which follow from the theories (Zhdanov and Trubnikov 
1989, Bychkov 1998a, Kazakov and Liberman 2002b, Joulin and Cambray 
1992). As we can see, the formulas (Zhdanov and Trubnikov 1989, Kazakov 
and Liberman 2002b, Joulin and Cambray 1992) overestimate the velocity 
increase noticeably, especially for 6 = 8 corresponding to stoichiometric 
methane and propane fiames. So far, the formula proposed in (Bychkov 
1998a) 

''ma. 2 93 + 62 + 36-1 ^ ^ 

provides the best analytical fit for the numerical results. At this point we 
have to note that, as was remarked in (Boury 2003), neither of the pa- 
pers (Zhdanov and Trubnikov 1989, Bychkov 1998a, Kazakov and Liberman 
2002b, Kazakov and Liberman 2002a) included Joulin counter-terms (see 
again (Joulin and Cambray 1992)), which would lead to considerable quan- 
titative corrections to all these results including equation (pTT|l . However, 
one can obtain almost any velocity increase in scope of the perturbation ap- 
proaches, and the formula (n?T|l of (Bychkov 1998a) as well as other formulas 
of (Zhdanov and Trubnikov 1989, Kazakov and Liberman 2002b, Kazakov 
and Liberman 2002a) may be equally treated as analytical guesses rather 
than unambiguous results. Taking (pTT| as the estimate for the velocity in- 
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crease we find 



The model equation involves also time-dependent nonlinear terms, which 
cannot be adjusted with the help of direct numerical simulations for station- 
ary flames. However, studies of curved flame stability (Bychkov, Kovalev and 
Liberman 1999, Petchenko and Bychkov 2004) indicate that time-dependent 
nonlinear terms are of minor importance. For simplicity, when constructing a 
qualitative model the time-dependent nonlinear term may be omitted. Still, 
there is another problem with formula (f3T|l . namely, equation (j32| does not 
reproduce correct asymptotics for Ci at 6 — 1. Making slight modifications 
of (n?T|l we can remedy this trouble, since, as we have pointed out above, the 
stationary flame velocity is almost a free parameter in the nonlinear pertur- 
bation theories. For example, we can choose 



_9 (9-1)^ 
^imax- 2 93 + 292 + 50-4 ^^^^ 



with respective corrections to the coefficient Ci 



The velocity increase (f33|) is shown in figure [T] by the dashed line. As we 
can see, it provides even better agreement with direct numerical simulations 
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than ((SH). 

Of course, the above calculations is not the way to construct an unam- 
biguous rigorous equation to calculate the flame propagation velocity. Par- 
ticularly, one of the referees points out that the second order time derivative 
is useless, if the ultimate goal is only to compute the flame velocity. But 
we recall that the non linear theory of the DL instability in an inevitable 
starting point for many other problems like oblique flames, tulip flames, 
flame interaction with turbulence, with acoustics or shock waves, burning 
in tubes with heat losses. As an example, the non linear equations (see for 
instance (Joulin and Cambray 1992, Bychkov 1998a)) developed to describe 
the DL instability were later used with some modifications to study turbu- 
lent flames in (Cambray and Joulin 1992, Denet 1997, Bychkov 2000, Zaytsev 
and Bychkov 2002, Akkerman and Bychkov 2003). In the same way, equa- 
tion ( (pHll ) can also be modifled to get an understanding of other, much more 
complicated phenomena. 

5 Conclusion 

In this article, starting from a low vorticity approach (Bychkov et al. 2003) 
proposed to describe premixed flames in a coordinate-free way, we have devel- 
oped this formulation in powers of the gas expansion parameter. It appears 
that at the lowest order in gas expansion, the Frankel equation (equivalent 
of the Sivashinsky equation for the coordinate-free case) is recovered. At the 
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second order, we have shown that complications arise in the case of oblique 
flames. On the contrary, if the tangential velocity is small (planar on aver- 
age and expanding flames), we have obtained a modified form of the Frankel 
equation, with a correction of the surface charge. In the planar case, we have 
shown that this modified equation reduces to the Sivashinsky-Clavin equa- 
tion (a second order in gas expansion correction to the original Sivashinsky 
equation). We have thus shown that the small vorticity formulation contains 
the Sivashinsky-Clavin equation as a particular case. A direct numerical so- 
lution of this formulation, although difficult, could describe both the slow 
dynamics of a fiame without external forcing, and the rapid evolution that 
takes place under some conditions (acoustic forcing, interaction with shock 
waves). On the contrary, the equations obtained here as the lowest orders of 
a development in gas expansion are inherently limited to a slow, first order 
in time dynamics. Although potentially easier to solve than the full small 
vorticity equations (at least without tangential blowing) it must be recalled 
that in the planar case, in order to obtain good quantitative agreement with 
numerical simulations, Joulin and Cambray (Joulin and Cambray 1992) have 
been obliged to perform some empirical modifications of the coefficients of the 
Sivashinsky-Clavin equation. The purpose of these new coefficients was to 
describe better the instability growth rates and the amplitude of stationary 
cellular fiames. We suggest different ways to construct similar phenomenolog- 
ical equations, particularly taking into account second-order time derivatives 
inherent to the DL dispersion relation. We also discuss the best way of ad- 
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justing the numerical coefficients of the model equation using recent results 
of direct numerical simulations for the velocity increase because of the DL 
instability (Bychkov et al. 1996, Kadowaki 1999, Travnikov et al. 2000). It 
remains to be seen if the same type of modification has to be used in the 
coordinate-free case, both for the small vorticity equations and its small gas 
expansion, low frequency limit. In any event, we hope that the present article 
has served to explain the relations between different existing approaches to 
the problem of nonlinear premixed flames dynamics. 
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Denet and Bychkov "Low vorticity and small gas expansion in premixed flames" 



Figure 1: Maximal velocity increase for curved stationary flames scaled by 
the planar flame velocity versus the thermal expansion B. The markers show 
results of direct numerical simulations (Bychkov et al. 1996, Travnikov et 
al. 2000) (circles) and (Kadowaki 1999) (crosses). The solid lines correspond 
to the analytical results of (Zhdanov and Trubnikov 1989, Bychkov 1998a, 
Kazakov and Liberman 2002b, Joulin and Cambray 1992). The dashed line 
presents equation 
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